!********************************************************************************************
!       THIS PROGRAM ACCOMPANIES Fried et al (2024)
!
!       "Understanding the Inequality and Welfare Impacts of Carbon Tax Policies"
!
!   This code is used to solve the steady state without the carbon tax
!   If you use this program, I would kindly ask that you cite my paper.
!*********************************************************************************************

PROGRAM hetero_ss
USE functions_mod
USE grid_interpolation
USE lambda2_optimizer_mod
USE retired_valuer_mod
USE working_valuer_mod
USE parameters_mod
USE parameters_energy_mod
USE NEQNF_INT

IMPLICIT NONE
INTERFACE
    function randu(iseed)
        INTEGER, intent(inout)   ::iseed
    DOUBLE PRECISION   ::randu
end function randu
END INTERFACE
include 'mpif.h'

parameter send_data_tag = 2001, return_data_tag = 2002

INTEGER            :: con_pos,kstart, capital_positive, tax_loop, non_clear,temp,kit,i,s,it1
INTEGER            :: ac,out_loop,it,e_it
INTEGER            :: sim_shocks(J,simulationsN),sim_types(J,simulationsN)
INTEGER            :: taxes_clear(max_outer_loops+1)
INTEGER            :: shocks_start1,shocks_start2
INTEGER            :: num_procs_run,seed,it2

INTEGER            :: my_id, root_process, ierr, status(MPI_STATUS_SIZE)
INTEGER            :: num_procs, an_id
INTEGER            :: num_procs_run2,optimal_baseline
INTEGER            :: start_row,end_row,increment,large_inc,increment2,large_inc2
INTEGER            :: bad_combination_inner_nostop,bad_combination_outer,bad_combination_inner_stop


DOUBLE PRECISION   :: b,pricefactor,C_for_share,Ce_for_share,xguess(3),solver_out(3),K1_share,N1_share,rate_use
DOUBLE PRECISION   :: curv,diff_parameter,diff_parameter_v(6),survive(J), dying(J),value,types(typesN)
DOUBLE PRECISION   :: growth,  pop_temp1(J),hours_jnr,sim_earnings(J,simulationsN),sim_capital_earnings(J,simulationsN)
DOUBLE PRECISION   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),sim_consumption(J,simulationsN)
DOUBLE PRECISION   :: sim_energy(J,simulationsN),sim_capital_tax(J,simulationsN),sim_labor_tax(J,simulationsN)
DOUBLE PRECISION   :: agg_labor(J,typesN),agg_capital(J+1,typesN),agg_consumption(J,typesN),agg_energy(J,typesN)
DOUBLE PRECISION   :: kprime(typesN,shocksN,kgridN,J-1),labor(typesN,shocksN,kgridN,Jnr),energy_ratio
DOUBLE PRECISION   :: dummy,energy_share_sim(J,simulationsN),sim_v(simulationsN)
DOUBLE PRECISION   :: v(typesN,shocksN,kgridN,J+1),v_retired(typesN,shocksN,kgridN,J+1),sim_consumption_tilda(J,simulationsN)
DOUBLE PRECISION   :: v_in(typesN,shocksN,kgridN),sim_utility(J,simulationsN),sim_consumption_dummy(1,simulationsN)
DOUBLE PRECISION   :: ctot80(kgridN,1),e80(kgridN,1),c80(kgridN,1), con(kgridN),aggregates(31),con_dummy
DOUBLE PRECISION   :: ss(typesN,(shocksN)/2),ss1(typesN,(shocksN)/2),aggregates_bigger(50)

DOUBLE PRECISION   :: energy_miss,N1_ss(typesN)
DOUBLE PRECISION   :: c_upper(kgridN),avg_earnings1,giniWelSS,lifeWelBar,giniadj,giniIncSS_aftertax
DOUBLE PRECISION   :: gs,increment_real,increment_real2

DOUBLE PRECISION   :: hbar,partprod, K_ig,E,Ec,Ec1,Ep, Ec_ig,N_ig, tr_ig, K, N, N1,N2,K1,K2,Ae,Atfp,sigmaa,gov_spend
DOUBLE PRECISION   :: K1_prime,N1_prime,Y, w, r,kgridl, kgridu,  tr1,lambda_ig,sigmav,sigmaeps,rho
DOUBLE PRECISION   :: a,a1,tax_diff,tax_diff_keep,govt_taxes,tol,pop_start
DOUBLE PRECISION   :: shocks_temp((shocksN)/2),stat_dist_temp((shocksN)/2),shock_pi_temp((shocksN)/2,(shocksN)/2)
DOUBLE PRECISION   :: shocks(shocksN),shock_pi(shocksN,shocksN,Jnr-1,typesN),shock_cumpi(shocksN,shocksN,J,typesN)
DOUBLE PRECISION   :: avg_earnings,avg_earnings_ig,prices(9)
DOUBLE PRECISION   :: avg_earnings_ss1(typesN,(shocksN)/2),avg_earnings_ss_ig(typesN,(shocksN)/2),avg_earnings_ss(typesN,(shocksN)/2)
DOUBLE PRECISION   :: random(J,simulationsN),skink(skinkN),sspay(skinkN)
DOUBLE PRECISION   :: ss_count(typesN,(shocksN)/2),sim_lifetime_after_tax_income(simulationsN)



DOUBLE PRECISION, allocatable, dimension(:,:) :: zeros, population,sim_labor_alloc,sim_capital_alloc,sim_capital_tax_alloc,sim_labor_tax_alloc
DOUBLE PRECISION, allocatable, dimension(:,:) :: sim_consumption_alloc,sim_earnings_alloc,sim_energy_alloc
DOUBLE PRECISION, allocatable, dimension(:) :: zeros2,sim_v_alloc
DOUBLE PRECISION, allocatable, dimension(:,:,:) :: v_alloc,kprime_alloc,labor_alloc


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MPI Parameters

root_process = 0
call MPI_INIT(ierr)
call MPI_COMM_RANK (MPI_COMM_WORLD, my_id, ierr)
call MPI_COMM_SIZE (MPI_COMM_WORLD, num_procs, ierr)
num_procs_run=num_procs
num_procs_run2=num_procs
!Increments for simulations
if (simulationsN<(num_procs-1)) then
    num_procs_run2=simulationsN+1
end if
increment_real2=simulationsN/(num_procs_run2-1.0D0)
increment2=floor(increment_real2)
large_inc2=anint((increment_real2-increment2)*(num_procs_run2-1))

!Increments for value function calculation

if (kgridN<(num_procs-1)) then
    num_procs_run=kgridN+1
end if
increment_real=kgridN/(num_procs_run-1.0D0)
increment=floor(increment_real)
large_inc=anint((increment_real-increment)*(num_procs_run-1))


CALL ERSET(4,1,0)
CALL ERSET(3,1,0)
CALL ERSET(2,1,0)

household=0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Computational Parameters
a=0.15D0
a1=0.15D0
tol=0.0005D0
allocate(zeros2(max_outer_loops+1))
zeros2=0
taxes_clear=zeros2
deallocate(zeros2)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Utility Parameters
sigma1=2.0D0


if(household==0) then
    !Individual
    gamma1=0.99074D0
    ebar=0.001324D0
!    ebar=0.00D0
    beta=0.99502D0
    chi=73.3
    sigma2=0.5D0
    kgridl=-0.157D0
elseif (household==1) then
    !Household
    gamma1=0.99275D0
    ebar=0.0068D0
    beta=1.0085D0
    chi=90.0D0
    sigma2=0.55D0
    kgridl=0.0D0
else
    print*,"Set Household or Individual"
end if

!!!1 if solving optimal, 0 if solving baseline
optimal_baseline=1

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Productivity Parameters
!Individual
hbar=0.00D0
partprod=1.0D0

if(household==0) then
epsilon1(:,1) = (/1.0D0, 1.05037955556458D0, 1.10168574952967D0, 1.15381030928342D0, 1.20663608757434D0, 1.26003734304819D0, 1.3138800980071D0, 1.36802257372016D0, 1.42231570281972D0, 1.47660371749329D0, 1.53072481133652D0, 1.584511871879D0, 1.63779327994293D0, 1.6903937711557D0, 1.74213535412359D0, 1.79283827899664D0, 1.84232204942572D0, 1.89040647024333D0, 1.93691272260084D0, 1.9816644577762D0, 2.02448890043734D0, 2.06521795181501D0, 2.10368928301267D0, 2.13974740856428D0, 2.1732447303485D0, 2.20404254208235D0, 2.23201198484835D0, 2.25703494445744D0, 2.27900488191109D0, 2.29782758879663D0, 2.31342186012357D0, 2.32572007787805D0, 2.33466869942757D0, 2.34022864584003D0, 2.34237558617614D0, 2.34110011486113D0, 2.33640782032644D0, 2.32831924422039D0, 2.31686973160448D0, 2.30210917366396D0, 2.28410164555347D0, 2.26292494305636D0, 2.23867002274592D0, 2.21144035128553D0, 2.18135117038086D0, 2.14852868468884D0/)
elseif(household==1) then
epsilon1(:,1) = (/1.0D0, 1.11972993442023D0, 1.23163188204661D0, 1.33628274802356D0, 1.43423225535758D0, 1.52600294491726D0, 1.61209017543323D0, 1.69296212349824D0, 1.76905978356708D0, 1.84079696795664D0, 1.90856030684589D0, 1.97270924827585D0, 2.03357605814965D0, 2.09146582023247D0, 2.14665643615159D0, 2.19939862539634D0, 2.24991592531814D0, 2.29840469113051D0, 2.34503409590901D0, 2.3899461305913D0, 2.4332556039771D0, 2.47505014272822D0, 2.51539019136855D0, 2.55430901228404D0, 2.59181268572273D0, 2.62788010979474D0, 2.66246300047226D0, 2.69548589158955D0, 2.72684613484296D0, 2.75641389979091D0, 2.7840321738539D0, 2.8095167623145D0, 2.83265628831737D0, 2.85321219286924D0, 2.87091873483891D0, 2.88548299095726D0, 2.89658485581726D0, 2.90387704187394D0, 2.90698507944442D0, 2.90550731670788D0, 2.89901491970559D0, 2.88705187234091D0, 2.86913497637924D0, 2.84475385144809D0, 2.81337093503703D0, 2.77442148249772D0/)
else
    print*,"Set Household"
end if
!Firm
alpha1 = 0.3D0
alpha2 = 1.0D0-0.403D0
theta=0.03D0
delta = 0.079062D0
Atfp=1.0D0
Ae=Atfp*1.0D0


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Heterogeneity Parameters
if(household==0) then
    rho=0.958D0
    sigmaa=0.065D0**0.5D0
    sigmav=(0.017D0/(1-rho**2))**0.5D0
    sigmaeps=0.081D0**0.5D0
elseif(household==1) then
    rho=0.972D0
    sigmaa=0.022D0**0.5D0
    sigmav=(0.012D0/(1-rho**2))**0.5D0
    sigmaeps=0.093D0**0.5D0
else
    print*,"Set Household"
end if



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Calculate Heterogeneity
do i =1,typesN
    types(i)=exp(-(sigmaa)+(real(i-1)/real(typesN-1))*2.0D0*sigmaa)
end do

! Stochastic parameters
CALL rouwenhorst(rho,sigmav,shock_pi_temp,shocks_temp,stat_dist_temp,(shocksN)/2)


shocks_temp=exp(shocks_temp)
shocks_temp=shocks_temp/(sum(stat_dist_temp*shocks_temp))

shocks(1:int(shocksN/2))=shocks_temp*exp(-sigmaeps)
shocks(int(shocksN/2)+1:int(shocksN))=shocks_temp*exp(sigmaeps)

shocks_start1=int((((shocksN)/2)+1)/2)
shocks_start2=int(shocks_start1+((shocksN)/2))

!E to E
do it1=1,Jnr-1
do it2=1,typesN
shock_pi(1:int((shocksN)/2),1:int((shocksN)/2),it1,it2)=shock_pi_temp/2.0D0
shock_pi(int((shocksN)/2+1):(shocksN),1:int((shocksN)/2),it1,it2)=shock_pi_temp/2.0D0
shock_pi(1:int((shocksN)/2),int((shocksN)/2+1):(shocksN),it1,it2)=shock_pi_temp/2.0D0
shock_pi(int((shocksN)/2+1):(shocksN),int((shocksN)/2+1):(shocksN),it1,it2)=shock_pi_temp/2.0D0
end do
end do






!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Population Parameters
growth = 0.011D0
if(household==0) then
pop_temp1 =(/197316.D0, 197141.D0, 196959.D0, 196770.D0, 196580.D0, 196392.D0, 196205.D0, 196019.D0, 195830.D0, 195634.D0, 195429.D0, 195211.D0, 194982.D0, 194739.D0, 194482.D0, 194211.D0, 193924.D0, 193619.D0, 193294.D0, 192945.D0, 192571.D0, 192169.D0, 191736.D0, 191271.D0, 190774.D0, 190243.D0, 189673.D0, 189060.D0, 188402.D0, 187699.D0, 186944.D0, 186133.D0, 185258.D0, 184313.D0, 183290.D0, 182181.D0, 180976.D0, 179665.D0, 178238.D0, 176689.D0, 175009.D0, 173187.D0, 171214.D0, 169064.D0, 166714.D0, 164147.D0, 161343.D0, 158304.D0, 155048.D0, 151604.D0, 147990.D0, 144189.D0, 140180.D0, 135960.D0, 131532.D0, 126888.D0, 122012.D0, 116888.D0, 111506.D0, 105861.D0, 99957.D0, 93806.D0, 87434.D0, 80882.D0, 74204.D0, 67462.D0, 60721.D0, 54053.D0,	47533.D0,	41241.D0, 35259.D0, 29663.D0, 24522.D0, 19890.D0, 15805.D0, 12284.D0, 9331.D0, 6924.D0, 5016.D0, 3550.D0, 2454.D0 /)
do i=1,(J-1)
   survive(i) = pop_temp1(i+1)/pop_temp1(i)
end do
survive(J) = 0.0D0
adult_equivalence=1.0D0
else
survive=(/1.0D0, 0.999443429520098D0, 0.999453048952094D0, 0.999474759951225D0, 0.999507517658172D0, 0.999550180294631D0, 0.99958688938571D0, 0.999616966794992D0, 0.999637259341625D0, 0.999646435563957D0, 0.999642137784891D0, 0.999638705313936D0, 0.999628264835928D0, 0.999611266702816D0, 0.999587623092204D0, 0.99955347402391D0, 0.999517484919676D0, 0.999474904683543D0, 0.999424216406889D0, 0.999368982231295D0, 0.999311637585569D0, 0.999247262201619D0, 0.99917475362204D0, 0.999096164795857D0, 0.999015366949988D0, 0.99892772761533D0, 0.998832359878012D0, 0.998729478390282D0, 0.998631176268071D0, 0.998535878125485D0, 0.998443264975948D0, 0.998338654631507D0, 0.998221906285934D0, 0.998078357864358D0, 0.997913954710353D0, 0.997725564043014D0, 0.997510605728885D0, 0.997277513805754D0, 0.99702940364917D0, 0.996758720419625D0, 0.996463887451456D0, 0.996136496024502D0, 0.995758248748826D0, 0.995343091935131D0, 0.994874730958277D0, 0.994354016314192D0, 0.993755337681475D0, 0.993065683870379D0, 0.992293047744931D0, 0.991433174915763D0, 0.990457448940283D0, 0.989303075213917D0, 0.987964810901374D0, 0.986475935412401D0, 0.984833077329247D0, 0.982973697702629D0, 0.980769120672249D0, 0.978166137633033D0, 0.975176536129091D0, 0.971754588528722D0, 0.967817497628447D0, 0.963253722024177D0, 0.957955278732916D0, 0.951933505633426D0, 0.945105096525417D0, 0.937431585348977D0, 0.92887640431547D0, 0.919444350437676D0, 0.909124748145856D0, 0.89794551325945D0, 0.885974030114678D0, 0.873304919941224D0, 0.859956433479889D0, 0.842665940409714D0, 0.824214308977887D0, 0.804792055183895D0, 0.785116734289799D0, 0.765756410717599D0, 0.747202819257837D0, 0.729458167443388D0, 0.713669733378468D0/)
survive(J) = 0.0D0
adult_equivalence=(/1.21998398833892D0, 1.26899136580154D0, 1.32056246723579D0, 1.3734138700101D0, 1.42641978111077D0, 1.47860215039009D0, 1.52912105574409D0, 1.57726536021961D0, 1.62244364105098D0, 1.66417539062616D0, 1.70208248938226D0, 1.73588095063074D0, 1.76537293731188D0, 1.79043905067891D0, 1.81103089091152D0, 1.82716388965886D0, 1.83891041451209D0, 1.84639314540638D0, 1.84977872295234D0, 1.84927166869702D0, 1.84510857731432D0, 1.83755258072496D0, 1.8268880841459D0, 1.81341577406915D0, 1.79744789817028D0, 1.77930381714618D0, 1.75930582848249D0, 1.73777526215036D0, 1.71502884823281D0, 1.69137535648057D0, 1.66711250779728D0, 1.64252415765435D0, 1.61787775143512D0, 1.5934220517087D0, 1.5693851374331D0, 1.54597267508807D0, 1.52336646173711D0, 1.50172324001931D0, 1.48117378507042D0, 1.46182226337349D0, 1.44374586353911D0, 1.42699469901482D0, 1.41159198272448D0, 1.39753447363654D0, 1.3847931952624D0, 1.37331442608375D0, 1.36302096190969D0, 1.35381365016324D0, 1.34557319609731D0, 1.33816224094031D0, 1.33142771197095D0, 1.32520344452275D0, 1.31931307591802D0, 1.31357321133145D0, 1.30779686158245D0, 1.30179715285813D0, 1.29539130836514D0, 1.28840490191053D0, 1.28067638341311D0, 1.27206187634362D0, 1.26244024709437D0, 1.25171844627878D0, 1.23983712196005D0, 1.22677650480938D0, 1.21256256519385D0, 1.19727344219358D0, 1.18104614454827D0, 1.16408352353373D0, 1.14666151776743D0, 1.12913666994348D0, 1.11195391549768D0, 1.09896503226485D0, 1.08988645284858D0, 1.0816144615019D0, 1.07430423674565D0, 1.0681188559846D0, 1.06322914430733D0, 1.0598135018135D0, 1.0580577094685D0, 1.05815471348597D0, 1.06030438823754D0/)
end if

dying = 1.D0 - survive
allocate(population(J,2))
population=0.0D0
population(1,1) = 1.0D0
do i=2,(J)
   population(i,1) = survive(i-1)*population((i-1),1)/(1+growth)
end do
pop_start=population(1,1)/sum(population(1:J,1))
population(:,2) = population(:,1)/sum(population(:,1))



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Grids

kgrid(1)=kgridl
kgrid(kgridnegN+1)=0.0D0
curv = 1.5D0
do ac=2,kgridnegN+1
	kgrid(kgridnegN-ac+2)=kgrid(kgridnegN+1)+(kgridl)*((ac-1.0D0)/(kgridnegN))**curv
end do

kgridu = 25.0D0
curv = 1.35D0
do ac=kgridnegN+2,kgridN
	kgrid(ac)=kgrid(kgridnegN+1)+(kgridu-kgrid(kgridnegN+1))*((ac-kgridnegN-1.0D0)/(kgridN-kgridnegN-1.0D0))**curv
end do
kgrid(kgridN)=kgridu
capital_positive = minloc(kgrid,1, mask = kgrid.ge.0)
kstart=capital_positive



!kgridl = 0.0D0
!kgridu = 101.0D0
!kstart = 1
!!Curved grids to handle outliers
!curv = 2.25D0
!kgrid(1)=kgridl
!
!do ac=2,kgridN
!	kgrid(ac)=kgrid(1)+(kgridu-kgridl)*((ac-1.0D0)/(kgridN-1.0D0))**curv
!end do
!kgrid(kgridN)=kgridu


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Social Security Payment Information
b=0.45D0
!!!!!!!!!!!!!!!!!!!!!!!!!!!! Government Parameters
gs=.157D0
govt_taxes=0.0D0
!lambda0=.258D0

if(optimal_baseline==0) then
!!!!Baseline calibration
lambda1=.031D0
lambdak=.36D0
else
!!!!Optimal policy
lambda1=.121D0
lambdak=.116D0
end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Initial Guesses
lambda_ig=0.826672024D0
K_ig=1.74964187009136D0
N_ig=0.525888026916112D0
Ec_ig=0.0084609597055167D0
!E_ig=16.5373251118981D0
tr_ig=0.0231357700461317D0
ss_ig(2,:)=(/0.224093152920768D0, 0.245707033732439D0, 0.271216751401476D0, 0.294198268700008D0, 0.305832866510048D0/)
ss_ig(1,:)=(/0.171185728929924D0, 0.187798546817684D0, 0.207066003615388D0, 0.229451373360864D0, 0.258223519766395D0/)
avg_earnings_ig=0.561330492585394D0
avg_earnings_ss_ig=0.644971442206701D0
avg_earnings_ss_ig(1,:)=(/0.328105717283857D0, 0.380020773183123D0, 0.440231575676112D0, 0.510185856130657D0, 0.60009881364819D0/)
avg_earnings_ss_ig(2,:)=(/0.493441417255325D0, 0.560984794791977D0, 0.640702662507457D0, 0.725536046370185D0, 0.803086275638889D0/)

!!!!!!!!!!!!!!!!!!!!!!!! Initialize loop stuff
tax_loop = 1
tax_diff=1
non_clear=0
bad_combination_inner_stop=0
bad_combination_inner_nostop=0
bad_combination_outer=0
diff_parameter=1.0D0


!Set Up population for simulations
seed=53415691
do s=1,simulationsN
do i=1,J
    random(i,s)=randu(seed)
end do
end do

i=1
do s=1,typesN
    sim_types(1,i:i+int(simulationsN/(typesN))-1)=s
if (my_id==1) then
print*,"i i+",i,i+int(simulationsN/(typesN))-1
end if
    i=i+int(simulationsN/(typesN))
end do

!sim_types(1,1:int((simulationsN)/2))=1
!sim_types(1,int((simulationsN)/2)+1:simulationsN)=2
do s=2,J
    sim_types(s,:)=sim_types(1,:)
end do

seed=52345691
do s=1,simulationsN
do i=1,J
    random(i,s)=randu(seed)
end do
end do

i=1
do it =1, typesN
    sim_shocks(1,i:i+int(simulationsN/(2*(typesN)))-1)=shocks_start1
if (my_id==1) then
print*,"i i+",i,i+int(simulationsN/(2*(typesN)))-1
end if
    i=i+int(simulationsN/(2*(typesN)))
    sim_shocks(1,i:i+int(simulationsN/(2*(typesN)))-1)=shocks_start2
if (my_id==1) then
print*,"i i+",i,i+int(simulationsN/(2*(typesN)))-1
end if
    i=i+int(simulationsN/(2*(typesN)))
end do

!
!sim_shocks(1,1:int((simulationsN)/4))=shocks_start1
!sim_shocks(1,int((simulationsN)/4)+1:int((simulationsN)/4)*2)=shocks_start2
!sim_shocks(1,int((simulationsN)/4)*2+1:int((simulationsN)/4)*3)=shocks_start1
!sim_shocks(1,int((simulationsN)/4)*3+1:simulationsN)=shocks_start2

do i=1,shocksN
do s=1,shocksN
do it1=1,Jnr-1
do it2=1,2
    shock_cumpi(i,s,it1,it2)=sum(shock_pi(i,1:s,it1,it2))
end do
end do
end do
end do

do it1=Jnr,J
do it2=1,2
    shock_cumpi(:,:,it1,it2)=shock_cumpi(:,:,Jnr-1,it2)
end do
end do

shock_cumpi(:,shocksN,:,:)=1.0D0


!Low
do s=1,int(simulationsN/2)
do i =2,J
    if (i<Jnr) then
        sim_shocks(i,s)=minloc(shock_cumpi(sim_shocks(i-1,s),:,i,1),1,mask=shock_cumpi(sim_shocks(i-1,s),:,i,1).ge.random(i,s))
    else
        sim_shocks(i,s)=sim_shocks(i-1,s)
    end if
end do
end do
!High
do s=int(simulationsN/2)+1,simulationsN
do i =2,J
    if (i<Jnr) then
        sim_shocks(i,s)=minloc(shock_cumpi(sim_shocks(i-1,s),:,i,2),1,mask=shock_cumpi(sim_shocks(i-1,s),:,i,2).ge.random(i,s))
    else
        sim_shocks(i,s)=sim_shocks(i-1,s)
    end if
end do
end do

if(my_id==1) then
print*,"shocks",shocks
end if


ss_count=0.0D0
do it =1 ,simulationsN
    if(sim_shocks(Jnr,it).le.shocksN/2) then
        ss_count(sim_types(1,it),sim_shocks(Jnr,it))=ss_count(sim_types(1,it),sim_shocks(Jnr,it))+1.0D0
    else
        ss_count(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)=ss_count(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)+1.0D0
    end if
end do
if(my_id==1) then
    print*,"ss count",ss_count(1,:)
    print*,"ss count",ss_count(2,:)
end if



if(optimal_baseline==0) then
if (my_id==1) then
    OPEN(UNIT=54, FILE="/<file path steady state no tax>/shocks.dat", ACTION="write", STATUS="replace", &
        FORM="formatted")
    WRITE(54,881) shocks
    close(unit=54)

    OPEN(UNIT=55, FILE="/<file path steady state no tax>/types.dat", ACTION="write", STATUS="replace", &
        FORM="formatted")
    WRITE(55,881) types
    close(unit=55)
end if
end if
881 format(F24.16)


out_loop=1
lambda0=lambda_ig
!capital_positive = minloc(kgrid,1, mask = kgrid.ge.0)

!MASTER
if (my_id .eq. root_process) then



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ENTER FUNCTION HERE IF MAKING A FUNCTION
!Initial guesses
K=K_ig
N=N_ig
tr=tr_ig
Ec=Ec_ig
!E=E_ig
ss=ss_ig
avg_earnings=avg_earnings_ig
avg_earnings_ss=avg_earnings_ss_ig
bad_combination_inner_stop=0
bad_combination_outer=0

!Loop for price convergence
DO WHILE (diff_parameter>tol)
   tax_loop = 1
   tax_diff=1.0D0
   non_clear=0

    passed_parameters(1)=K
    passed_parameters(2)=N
    passed_parameters(3)=Ec
    passed_parameters(4)=Atfp
    passed_parameters(5)=Ae
    passed_parameters(6)=0.0D0

    if(out_loop==1) then
        xguess(1)=0.05D0
    else
        xguess(1)=Ep
    end if
    if(out_loop==1) then
        xguess(2)=K*.5D0
        xguess(3)=N*.6D0
    else
        xguess(2)=K*K1_share
        xguess(3)=N*N1_share
    end if

    print*,"parms",Atfp,Ae,xguess,Ec
    CALL D_NEQNF(energy_root_finder,solver_out,ERRREL=0.0001D0,xguess=xguess)
    Ep=solver_out(1)
    E=Ep+Ec
    K1=solver_out(2)
    K2=K-K1
    N1=solver_out(3)
    N2=N-N1
    K1_share=K1/K
    N1_share=N1/N

    Y=Atfp*(K1**(alpha1)*N1**(1.0D0-alpha1-theta))*Ep**theta
    pe=Atfp*theta*(K1/Ep)**alpha1*(N1/Ep)**(1.0D0-alpha1-theta)
    gov_spend=(Y+pe*Ec)*gs
    w=Atfp*(1.0D0-alpha1-theta)*(K1/N1)**(alpha1)*(Ep/N1)**theta
    r=Atfp*alpha1*(N1/K1)**(1.0D0-alpha1-theta)*(Ep/K1)**theta-delta
    tauss=sum(sum(ss*ss_count,2),1)/real(simulationsN)*sum(population(Jnr:J,2))/&
        (sum(sum(avg_earnings_ss*ss_count,2),1)/real(simulationsN)*sum(population(1:Jnr-1,2)))
!print*,"tauss",sum(ss)/real(typesN)*sum(population(Jnr:J,2))/(w*N),sum(ss)/real(typesN),sum(population(Jnr:J,2)),(w*N),ss(1),ss(2)
    avg_earnings_scale=avg_earnings*(1.0D0-tauss/2.0D0)
    ss_cap=avg_earnings*2.42

print*,"prices",w,r,avg_earnings_scale
!Extra loop for convergence on tax parameters (if turned on need tax parameter convergence begore prices update)
    DO WHILE (tax_diff >tol)
!Solve for policy funcation
        v = -1.0D0/0.0D0
        v(:,:,:,J+1)=0.0D0
        v_retired = -1.0D0/0.0D0
        kprime=0
        v_retired(:,:,capital_positive:kgridN,J+1)=0.0D0

!Last generation
            do it=capital_positive,kgridN
                do it2=1,typesN
                    do e_it=1,shocksN
                        if(e_it .le. (shocksN)/2) then
                            ctot80(it,1)=retcon(kgrid(it),0.0D0,tr,r,ss(it2,e_it),avg_earnings)
                            energy_miss=energy_sharer(0.00D0,ctot80(it,1)*.15D0,ctot80(it,1),energy_share_equation,.0001D0,c80(it,1),ctot80(it,1),J)
                            e80(it,1)=(ctot80(it,1)-c80(it,1))/pe
                            v_retired(it2,e_it,it,J) =utility_last(c80(it,1),e80(it,1),adult_equivalence(J))
                        else
                            v_retired(it2,e_it,it,J) =v_retired(it2,e_it-(shocksN/2),it,J)
                        end if
                    end do
                end do
            end do
!print*,"v_retired",J,v_retired(1,1:5,31:33,J)
!print*,"v_retired",s,v_retired(2,:,20,J),capital_positive
!print*,"v_retired",J,v_retired(1,1,20,J)

!Retired generations
        allocate(zeros(kgridN,1))
        zeros=0.0D0
        do s =(J-1),Jnr, -1
            do e_it=1,shocksN
            if(e_it .le. (shocksN)/2) then
                do it =1,typesN
                do kit=capital_positive, kgridN
                    con=-1.0D0
                    con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,kgrid(kit),kgrid(capital_positive:),&
                        tr,r,ss(it,e_it),avg_earnings)
                    con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))
                    if (con(kgridN)>0.015+ebar) then
                        con_pos=kgridN
                    end if
                    if (con(capital_positive+1)<=0.015D0+ebar) then
                        kprime(it,e_it,kit,s)=0
                        v_retired(it,e_it,kit,s)=-1.0D0/0.0D0
                    else
                    CALL retired_valuer(kgrid(kit),0.0D0,kgrid(con_pos),ss(it,e_it),tr,r,survive(s),v_retired(it,e_it,:,s+1),&
                        v_retired(it,e_it,kit,s),kprime(it,e_it,kit,s),avg_earnings,s)
                    end if
               end do
               end do
            else
                do it=1,typesN
                do kit=capital_positive, kgridN
                    kprime(it,e_it,kit,s)=kprime(it,e_it-(shocksN/2),kit,s)
                    v_retired(it,e_it,kit,s)=v_retired(it,e_it-(shocksN/2),kit,s)
                end do
                end do
            end if
            end do
!print*,"v_retired",s,v_retired(2,:,20,s)
!print*,"v_retired",s,v_retired(1,1,20,s)
!
        end do
        deallocate(zeros)
!        do e_it=1,shocksN
!            do i=1,typesN
!                v_retired(i,e_it,:,Jnr:J-1)=v_retired(1,1,:,Jnr:J-1)
!                    kprime(i,e_it,:,Jnr:J-1)=kprime(1,1,:,Jnr:J-1)
!            end do
!        end do



!Working generations
v(:,:,:,Jnr)=v_retired(:,:,:,Jnr)
                prices(1)=tr
                prices(2)=lambda0
                prices(3)=w
                prices(4)=r
                prices(5)=tauss
!                prices(6)=ss
                prices(7)=ss_cap
                prices(8)=avg_earnings
                prices(9)=pe

                CALL MPI_Bcast( prices, 9, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
                CALL MPI_Bcast( ss, typesN*(shocksN/2), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
        do s =(Jnr-1),2,-1
!print*,"workings",s
            CALL MPI_Bcast( v(:,:,:,s+1),typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            end_row=0
            do an_id = 1, num_procs_run -1
                start_row=end_row+1
                if (an_id>large_inc) then
                    end_row=start_row+increment-1
                else
                    end_row=start_row+increment
                end if
                call MPI_SEND(start_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
                call MPI_SEND(end_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
            end do
            do an_id=1,num_procs_run -1
                call MPI_RECV( start_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( end_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( labor(:,:,start_row:end_row,s),typesN*shocksN*(end_row-start_row+1) , MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( kprime(:,:,start_row:end_row,s),typesN*shocksN*(end_row-start_row+1) , MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( v(:,:,start_row:end_row,s),typesN*shocksN*(end_row-start_row+1) , MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
            end do
!print*,"v",s,v(2,:,40,s)
!print*,"kprime",s,kprime(2,:,40,s)
!print*,"labor",s,labor(2,:,40,s)
    !print*,"v working",s,v(2,3,20,s)
        end do




!! Simulate peoples

    CALL MPI_Bcast( v(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    CALL MPI_Bcast( v_retired(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    end_row=0
    do an_id = 1, num_procs_run2 -1
        start_row=end_row+1
        if (an_id>large_inc2) then
            end_row=start_row+increment2-1
        else
            end_row=start_row+increment2
        end if
        call MPI_SEND(start_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
        call MPI_SEND(end_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
    end do
    do an_id=1,num_procs_run2 -1
        call MPI_RECV( start_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( end_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_labor(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_earnings(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_capital(:,start_row:end_row),(J+1)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_consumption(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_energy(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_labor_tax(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_capital_tax(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_v(start_row:end_row),(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
    end do


print*,"value",sum(sim_v)/real(simulationsN)
print*,"max capital",maxval(sim_capital)
!print*,"lab",sim_labor(:,1)
!print*,"ear",sim_earnings(:,1)
!print*,"con",sim_consumption(:,1)
!print*,"eng",sim_energy(:,1)
!print*,"cap",sim_capital(:,1)





!Calculate Aggregate Profiles

agg_labor(:,1)=sum(sim_labor(:,1:int((simulationsN)/2)),2)/int((simulationsN)/2.0D0)
agg_labor(:,2)=sum(sim_labor(:,int((simulationsN)/2)+1:simulationsN),2)/(simulationsN-int((simulationsN)/2.0D0))
agg_capital(:,1)=sum(sim_capital(:,1:int((simulationsN)/2)),2)/int((simulationsN)/2.0D0)
agg_capital(:,2)=sum(sim_capital(:,int((simulationsN)/2)+1:simulationsN),2)/(simulationsN-int((simulationsN)/2.0D0))
agg_consumption(:,1)=sum(sim_consumption(:,1:int((simulationsN)/2)),2)/int((simulationsN)/2.0D0)
agg_consumption(:,2)=sum(sim_consumption(:,int((simulationsN)/2)+1:simulationsN),2)/(simulationsN-int((simulationsN)/2.0D0))
agg_energy(:,1)=sum(sim_energy(:,1:int((simulationsN)/2)),2)/int((simulationsN)/2.0D0)
agg_energy(:,2)=sum(sim_energy(:,int((simulationsN)/2)+1:simulationsN),2)/(simulationsN-int((simulationsN)/2.0D0))
!print*,"sim consumption1",agg_consumption(:,1)
!print*,"sim consumption2",agg_consumption(:,2)
!print*,"sim labor1",agg_labor(:,1)
!print*,"sim labor2",agg_labor(:,2)
!print*,"sim capital1",agg_capital(:,1)
!print*,"sim capital2",agg_capital(:,2)
!
!print*,"cap",agg_capital(:,1)
!print*,"con",agg_consumption(:,1)
!print*,"labr",agg_labor(:,1)
!
!print*,"scap",sim_capital(:,1)
!print*,"scon",sim_consumption(:,1)
!print*,"slabr",sim_labor(:,1)


!!Find lambda 2 that clears the market
CALL lambda2_optimizer(lambda0_1,gov_spend,sim_labor,sim_shocks,sim_types,sim_capital,hbar,w,shocks,&
                    types,partprod,r,tr,population(:,2),avg_earnings,ss,survive)

        tax_diff=abs((lambda0_1)-lambda0)/lambda0
        tax_diff_keep=tax_diff
        lambda0=(a1*lambda0_1+(1-a1)*lambda0)
        govt_taxes=tax_differ2(lambda0,sim_labor,sim_shocks,sim_types,sim_capital,hbar,w,shocks,types,partprod,&
                    r,tr,0.0D0,population(:,2),avg_earnings,ss,survive)
print*,"govt_taxes",govt_taxes,gov_spend
        govt_taxes=tax_differ2(lambda0,sim_labor,sim_shocks,sim_types,sim_capital,hbar,w,shocks,types,partprod,&
                    r,tr,gov_spend,population(:,2),avg_earnings,ss,survive)

        !Check for bad convergence within tax loop
!!        if (value==-1000000) then
!!            tax_diff=0
!!            bad_combination_inner_stop=1
!!        end if

        tax_loop=tax_loop+1
        if (tax_loop > max_tax_loop) then
          tax_diff = 0
          bad_combination_inner_nostop=1
        end if
        CALL MPI_Bcast( tax_diff, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    END DO




    K1_prime=0.0D0
    N1_prime=0.0D0
    N1_ss=0.0D0
    tr1=0.0D0
    hours_jnr=0.0D0
    Ec1=0.0D0

    C_for_share=0.0D0
    Ce_for_share=0.0D0
    avg_earnings1=0.0D0
    energy_ratio=0.0D0
    energy_share_sim=0.0D0
    avg_earnings_ss1=0.0D0
    do it=1,simulationsN
        do s=1,J
            C_for_share=C_for_share+sim_consumption(s,it)*population(s,2)/real(simulationsN)
            Ce_for_share=Ce_for_share+sim_energy(s,it)*population(s,2)/real(simulationsN)
            energy_share_sim(s,it)=sim_energy(s,it)*pe/((sim_energy(s,it)*pe)+sim_consumption(s,it))*population(s,2)
!            energy_ratio=energy_ratio+energy_share_sim(s,it)/real(simulationsN)
        end do
    end do
energy_ratio=Ce_for_share*pe/(C_for_share+Ce_for_share*pe)

    !Calculate new aggregates
    do it=1,simulationsN
        do s=1,Jnr-1
            avg_earnings1=avg_earnings1+sim_earnings(s,it)*population(s,2)/(simulationsN*sum(population(1:Jnr-1,2)))
            if (sim_shocks(Jnr,it) .le. shocksN/2) then
                avg_earnings_ss1(sim_types(1,it),sim_shocks(Jnr,it))=&
                    avg_earnings_ss1(sim_types(1,it),sim_shocks(Jnr,it))+&
                    min(ss_cap,sim_earnings(s,it))*population(s,2)/(ss_count(sim_types(1,it),sim_shocks(Jnr,it))*&
                    sum(population(1:Jnr-1,2)))
            else
                avg_earnings_ss1(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)=&
                    avg_earnings_ss1(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)+&
                    min(ss_cap,sim_earnings(s,it))*population(s,2)/(ss_count(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)*&
                    sum(population(1:Jnr-1,2)))
            end if
            sim_capital_earnings(s,it)=(sim_capital(s,it)+tr)*r
            Ec1=Ec1+sim_energy(s,it)*population(s,2)/(simulationsN)
        end do
        do s=1,Jnr-1
            if (sim_labor(s,it)<hbar) then
                N1_prime=N1_prime+partprod*sim_labor(s,it)*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*population(s,2)/simulationsN
                N1_ss(sim_types(s,it))=N1_ss(sim_types(s,it))+partprod*sim_labor(s,it)*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*population(s,2)/simulationsN*real(typesN)
            else
                N1_prime=N1_prime+sim_labor(s,it)*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*population(s,2)/simulationsN
                N1_ss(sim_types(s,it))=N1_ss(sim_types(s,it))+sim_labor(s,it)*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*population(s,2)/simulationsN*real(typesN)
            end if
            K1_prime=K1_prime+sim_capital(s+1,it)*population(s,2)/simulationsN
            tr1=tr1+sim_capital(s+1,it)*population(s,2)/simulationsN*(1.0D0-survive(s))
            hours_jnr=hours_jnr+sim_labor(s,it)*population(s,2)/simulationsN
        end do
        do s=Jnr,J
            Ec1=Ec1+sim_energy(s,it)*population(s,2)/(simulationsN)
            sim_capital_earnings(s,it)=(sim_capital(s,it)+tr)*r
            K1_prime=K1_prime+sim_capital(s+1,it)*population(s,2)/simulationsN
            tr1=tr1+sim_capital(s+1,it)*population(s,2)/simulationsN*(1.0D0-survive(s))
        end do
    end do


!        do it=1,typesN
!        end do
        hours_jnr=hours_jnr/sum(population(1:Jnr-1,2))
        K1_prime=K1_prime/(1.0D0+growth)
        tr1=tr1/(1.0D0+growth)

print*,"capital check",sum(sum(sim_capital(1:J,:),2)/real(simulationsN)*population(1:J,2))+tr1,K1_prime
    diff_parameter_v(1)=abs(((tr1)-tr))/tr
    diff_parameter_v(2)=abs(((K1_prime)-K))/K
    diff_parameter_v(3)=abs(((N1_prime)-N))/N
!    diff_parameter_v(4)=abs(((ss1*a+ss*(1-a))-ss))/ss
!    abs(((ss1(2))-ss(2)))/ss(2))
    diff_parameter_v(5)=abs(tax_diff_keep)
    print*,"K1",K1_prime,K
    print*,"N1",N1_prime,N
    print*,"tr1",tr1,tr
    print*,"ss1",ss1,ss
    print*,"hours",hours_jnr
!!!!Update Aggregates
    Ec=(Ec1*a+Ec*(1-a))
    N=(N1_prime*a+N*(1-a))
    K=(K1_prime*a+K*(1-a))


    tr=(tr1*a+tr*(1-a))
    avg_earnings=(avg_earnings1*a+avg_earnings*(1-a))
    avg_earnings_ss=(avg_earnings_ss1*a+avg_earnings_ss*(1-a))
    print*,"avg earnings ss",avg_earnings_ss

    skink(1)=0.0D0
    skink(2)=0.21D0*sum(avg_earnings_ss*ss_count)/sum(ss_count)
    skink(3)=1.29D0*sum(avg_earnings_ss*ss_count)/sum(ss_count)
    skink(4)=2.42D0*sum(avg_earnings_ss*ss_count)/sum(ss_count)
    skink(5)=10000000.0D0*1.4D0

    sspay(1)=0.0D0
    sspay(2)=skink(2)*0.90D0
    sspay(3)=sspay(2)+(skink(3)-skink(2))*0.32D0
    sspay(4)=sspay(3)+(skink(4)-skink(3))*0.15D0
    sspay(5)=sspay(4)
!    do it=1,sgridN
!        spayoutgrid(it)=spayer(skinkN,skink,sspay,sgrid(it))*retire_scaler(retire_age)
!    end do
    do e_it=1,shocksN/2
        do it=1,typesN
            ss1(it,e_it)=spayer(skinkN,skink,sspay,avg_earnings_ss1(it,e_it))
        end do
    end do
    diff_parameter_v(4)=maxval(abs(((ss1)-ss))/ss)
    diff_parameter_v(6)=max(abs(((avg_earnings1)-avg_earnings))/avg_earnings,abs(((Ec1)-Ec))/Ec,&
        maxval(abs(((avg_earnings_ss1)-avg_earnings_ss))/avg_earnings_ss1))
    diff_parameter=MAXVAL(diff_parameter_v)

    print*,"ss1",ss1(1,:),ss1(2,:)
    print*,"ss",avg_earnings_ss1(1,:),avg_earnings_ss1(2,:),sum(avg_earnings_ss1*ss_count)/sum(ss_count)

    do e_it=1,shocksN/2
        do it=1,typesN
            ss(it,e_it)=(ss1(it,e_it)*a+ss(it,e_it)*(1-a))
        end do
    end do

    if (bad_combination_inner_stop == 1) then
        diff_parameter=0.0D0
    end if
    taxes_clear(out_loop)=0
    if (abs(govt_taxes/gov_spend)>tol) then
        taxes_clear(out_loop)=1
        non_clear=1
        diff_parameter=1.0D0
    end if
    if (taxes_clear(out_loop)==0) then
        allocate(zeros2(max_outer_loops+1))
        zeros2=0
        taxes_clear=zeros2
        deallocate(zeros2)
    end if
    if (sum(taxes_clear)>max_non_clear) then
        diff_parameter=0.0D0
        bad_combination_inner_stop=1
    end if
    if (out_loop > max_outer_loops) then
        bad_combination_outer=1
        diff_parameter=0.0D0
    end if
    out_loop=out_loop+1
    print*,"Loop"
    print*,"Plus one"
    print*,out_loop

    print*,"lambda0",lambda0
    print*,"govt miss",govt_taxes/gov_spend
    print*,"diffparm",diff_parameter
    print*,"diffvec",diff_parameter_v
    print*,"taxdiff",tax_diff_keep
    print*,"E",E
    print*,"Ec",Ec
    print*,"K",K
    print*,"N",N
    print*,"Y",Y
    print*,"ss",ss(1,:)
    print*,"ss",ss(2,:)
    print*,"tr",tr
    print*,"lambda",lambda0
    print*,"tr",tr
    print*,"Y",Y
    print*,"C",C_for_share
    print*,"K",K
    print*,"gov",gov_spend
    print*,"resource const",Y-C_for_share-gov_spend-K*delta
    print*,"resource const",Y-C_for_share-gov_spend-K*(delta+growth)
    print*,"r",r
    print*,"net r",r*(1-lambdak)
    print*,"max capital",maxval(sim_capital)
    print*,"avg_earnings",avg_earnings
    print*,"avg_earnings_ss",avg_earnings_ss
    print*,bad_combination_inner_nostop
    CALL MPI_Bcast( diff_parameter, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )

!print*,"con tilda",sum(sim_consumption_tilda,2)/real(simulationsN)
!print*,"con connn",sum(sim_consumption,2)/real(simulationsN)
!print*,"con engen",sum(sim_energy,2)/real(simulationsN)
!print*,"con labor",sum(sim_labor,2)/real(simulationsN)
!print*,"con savings",sum(sim_capital,2)/real(simulationsN)



END DO

do it =J,1,-1
    do it2=1,simulationsN
        if (it==J) then
            sim_utility(it,it2)=utility_last(sim_consumption(it,it2),sim_energy(it,it2),adult_equivalence(it))
        elseif(it>Jnr-1) then
            sim_utility(it,it2)=valuer_ret(sim_consumption(it,it2),sim_energy(it,it2),survive(it),sim_utility(it+1,it2),adult_equivalence(it))
        else
            sim_utility(it,it2)=valuer(sim_consumption(it,it2),sim_energy(it,it2),sim_labor(it,it2),survive(it),sim_utility(it+1,it2),adult_equivalence(it))
        end if
            sim_consumption_tilda(it,it2)=consumption_tilda(sim_consumption(it,it2),sim_energy(it,it2),adult_equivalence(it))
    end do
end do


sim_lifetime_after_tax_income=0.0D0
do it =J,1,-1
    do it2=1,simulationsN
        if (it==J) then
            sim_lifetime_after_tax_income(it2)=sim_earnings(it,it2)+sim_capital_earnings(it,it2)-sim_capital_tax(it,it2)-sim_labor_tax(it,it2)
        elseif(it>Jnr-1) then
            sim_lifetime_after_tax_income(it2)=1.0D0/(1.0D0+r*(1-lambdak))*sim_lifetime_after_tax_income(it2)+sim_earnings(it,it2)+sim_capital_earnings(it,it2)-sim_capital_tax(it,it2)-sim_labor_tax(it,it2)
        else
            sim_lifetime_after_tax_income(it2)=1.0D0/(1.0D0+r*(1-lambdak))*sim_lifetime_after_tax_income(it2)+sim_earnings(it,it2)+sim_capital_earnings(it,it2)-sim_capital_tax(it,it2)-sim_labor_tax(it,it2)
        end if


    end do
end do

!Gini for after tax income
giniIncSS_aftertax=0.0D0
lifeWelBar=sum(sim_lifetime_after_tax_income(:))/real(simulationsN)
giniadj=abs(2.0D0*real(simulationsN)**2.0D0*lifeWelBar)
do it =1,simulationsN
    do it2=1,simulationsN
        giniIncSS_aftertax=giniIncSS_aftertax+abs(sim_lifetime_after_tax_income(it)-sim_lifetime_after_tax_income(it2))/giniadj
    end do
end do


giniWelSS=0.0D0
lifeWelBar=sum(sim_utility(1,:))/real(simulationsN)
giniadj=abs(2.0D0*real(simulationsN)**2.0D0*lifeWelBar)

do it =1,simulationsN
    do it2=1,simulationsN
        giniWelSS=giniWelSS+abs(sim_utility(1,it)-sim_utility(1,it2))/giniadj
    end do
end do

print*,"gini",giniWelSS


do it =J,1,-1
    do it2=1,simulationsN
        if (it==J) then
            sim_consumption_dummy(1,it2)=sim_consumption(it,it2)


        else
            sim_consumption_dummy(1,it2)=sim_consumption(it,it2)+survive(it)*beta*sim_consumption_dummy(1,it2)
        end if
    end do
end do



aggregates(1)=w
aggregates(2)=r
aggregates(3)=K
aggregates(4)=N
aggregates(5)=tr
aggregates(6)=lambda0
aggregates(7)=gov_spend
aggregates(8)=tauss
aggregates(9)=maxval(sim_capital)
aggregates(10)=maxval(sim_capital)
aggregates(11)=E
aggregates(12)=lambdak
aggregates(13)=lambda1
aggregates(14)=Ce_for_share
aggregates(15)=C_for_share
aggregates(16)=pe
aggregates(17)=avg_earnings
aggregates(18)=giniWelSS
aggregates(19)=sum(sim_v)/real(simulationsN)
aggregates(20)=sum(sim_consumption_dummy(1,:))/real(simulationsN)
aggregates(21)=hours_jnr
aggregates(22)=Ec
aggregates(23)=K1
aggregates(24)=N1
aggregates(25)=K2
aggregates(26)=N2
aggregates(27)=Ep
aggregates(28)=maxval(sim_capital)
aggregates(29)=sum(sum(sim_consumption(:,:),2)*population(1:J,2))
aggregates(30)=sum(sum(sim_consumption_tilda(:,:),2)*population(1:J,2))
aggregates(31)=sum(sum(sim_energy(:,:),2)*population(1:J,2))

aggregates_bigger=0.0D0
aggregates_bigger(1:31)=aggregates(1:31)
aggregates_bigger(32)=giniIncSS_aftertax
aggregates_bigger(33)=alpha1
aggregates_bigger(34)=alpha2
aggregates_bigger(35)=theta
aggregates_bigger(36)=delta
aggregates_bigger(37)=Atfp
aggregates_bigger(38)=Ae
aggregates_bigger(39)=gamma1
aggregates_bigger(40)=ebar
aggregates_bigger(41)=beta
aggregates_bigger(42)=chi
aggregates_bigger(43)=sigma2
aggregates_bigger(44)=kgridl
aggregates_bigger(45)=lambda1
aggregates_bigger(46)=lambdak



         OPEN(UNIT=14, FILE="/<file name output>/v_retired_output.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(14,887) v_retired
          CLOSE(UNIT=14)


         OPEN(UNIT=15, FILE="/<file name output>/sim_labor_output.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(15,887) sim_labor
          CLOSE(UNIT=15)

          OPEN(UNIT=16, FILE="/<file name output>/sim_energy_output.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(16,887) sim_energy
          CLOSE(UNIT=16)

         OPEN(UNIT=17, FILE="/<file name output>/sim_capital_output.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(17,887) sim_capital
          CLOSE(UNIT=17)

         OPEN(UNIT=18, FILE="/<file name output>/sim_consumption_output.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) sim_consumption
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<file name output>/sim_utility_output.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) sim_utility
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<file name output>/sim_capital_tax_output.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) sim_capital_tax
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<file name output>/sim_labor_tax_output.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) sim_labor_tax
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<file name output>/sim_labor_earnings_output.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) sim_earnings
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<file name output>/sim_capital_earnings_output.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) sim_capital_earnings
          CLOSE(UNIT=18)

         OPEN(UNIT=13, FILE="/<file name output>/v_output.dat",STATUS="replace", &
               FORM="formatted")
          WRITE(13,887) v
          CLOSE(UNIT=13)

         OPEN(UNIT=19, FILE="/<file name output>/aggregates_output.dat",STATUS="replace", &
               FORM="formatted")
          WRITE(19,887) aggregates
          CLOSE(UNIT=13)

         OPEN(UNIT=19, FILE="/<file name output>/aggregates_bigger.dat",STATUS="replace", &
               FORM="formatted")
          WRITE(19,887) aggregates_bigger
          CLOSE(UNIT=13)

         OPEN(UNIT=19, FILE="/<file name output>/ss_output.dat",STATUS="replace", &
               FORM="formatted")
          WRITE(19,887) ss
          CLOSE(UNIT=13)

         OPEN(UNIT=19, FILE="/<file name output>/avg_earnings_ss_output.dat",STATUS="replace", &
               FORM="formatted")
          WRITE(19,887) avg_earnings_ss
          CLOSE(UNIT=13)



        OPEN(UNIT=21, FILE="/<file name output>/sim_shocks_output.dat",STATUS="replace", &
               FORM="formatted")
          WRITE(21,886) sim_shocks
          CLOSE(UNIT=21)

887 format(F24.16)
886 format(i9)

889 format(2F24.16)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Slaves
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!Slaves that are going to work
else if (my_id < (num_procs_run)) then
DO WHILE (diff_parameter>tol)
tax_diff=1.0D0
    DO WHILE (tax_diff >tol)
            CALL MPI_Bcast( prices, 9, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            CALL MPI_Bcast( ss, typesN*(shocksN/2), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            tr=prices(1)
            lambda0=prices(2)
            w=prices(3)
            r=prices(4)
            tauss=prices(5)
!            ss=prices(6)
            ss_cap=prices(7)
            avg_earnings=prices(8)
            pe=prices(9)
            avg_earnings_scale=avg_earnings*(1.0D0-tauss/2.0D0)

!Solve Working
        do s=Jnr-1,2,-1
            CALL MPI_Bcast( v_in(:,:,:), typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            call MPI_RECV(start_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
            call MPI_RECV( end_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
            allocate(labor_alloc(typesN,shocksN,end_row-start_row+1))
            allocate(kprime_alloc(typesN,shocksN,end_row-start_row+1))
            allocate(v_alloc(typesN,shocksN,end_row-start_row+1))
            labor_alloc=0.0D0
            kprime_alloc=0
            v_alloc=0.0D0
            do kit=start_row,end_row
                if(kgrid(kit)<0.0D0) then
                    rate_use=r/survive(s-1)
                else
                    rate_use=r
                end if
                do e_it=1,shocksN
                    do i=1,typesN
                        c_upper=c_upper_fnc(epsilon1(s,1),kgrid(kit),kgrid,kgridN,tr,&
                            w*shocks(e_it)*types(i),rate_use)
                        temp=int(MINLOC(c_upper,1,mask=c_upper.gt.0.015D0+ebar))
                        if (c_upper(kgridN)>0.015D0+ebar) then
                            temp = kgridN
                        end if
                        if (c_upper(2)<=0.015D0+ebar) then
                            kprime_alloc(i,e_it,kit-start_row+1)=0
                            v_alloc(i,e_it,kit-start_row+1)=-1.0D0/0.0D0
                        else
                            CALL working_valuer(s,kgrid(kit),temp,tr,rate_use,&
                            survive(s),&
                            w*shocks(e_it)*types(i)*epsilon1(s,1),shock_pi(e_it,:,s,i),&
                            v_in(i,:,1:temp),&
                            v_alloc(i,e_it,kit-start_row+1),kprime_alloc(i,e_it,kit-start_row+1),&
                            labor_alloc(i,e_it,kit-start_row+1),&
                            partprod,hbar)
                        end if
                    end do
                end do
            end do
            call MPI_SEND( start_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( end_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( labor_alloc,typesN*shocksN*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( kprime_alloc,typesN*shocksN*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( v_alloc,typesN*shocksN*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
            deallocate(labor_alloc)
            deallocate(kprime_alloc)
            deallocate(v_alloc)
        end do


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Simulating
    CALL MPI_Bcast( v(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    CALL MPI_Bcast( v_retired(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    call MPI_RECV(start_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
    call MPI_RECV(end_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

    allocate(sim_labor_alloc(J,end_row-start_row+1))
    allocate(sim_earnings_alloc(J,end_row-start_row+1))
    allocate(sim_capital_alloc(J+1,end_row-start_row+1))
    allocate(sim_consumption_alloc(J,end_row-start_row+1))
    allocate(sim_energy_alloc(J,end_row-start_row+1))
    allocate(sim_capital_tax_alloc(J,end_row-start_row+1))
    allocate(sim_labor_tax_alloc(J,end_row-start_row+1))
    allocate(sim_v_alloc(end_row-start_row+1))



    !!Initialize

    sim_capital_alloc(1,:)=0.0D0
    sim_labor_alloc=0.0D0
    sim_earnings_alloc=0.0D0
    sim_capital_tax_alloc=0.0D0
    sim_labor_tax_alloc=0.0D0





    do it=start_row,end_row
        !Working
        do s=1,Jnr-1
            if(sim_capital_alloc(s,it-start_row+1)<0.0D0) then
                if(s==1) then
                    rate_use=r/survive(1)
                else
                    rate_use=r/survive(s-1)
                end if
            else
                rate_use=r
            end if
            CALL working_valuer(s,sim_capital_alloc(s,it-start_row+1),kgridN,tr,rate_use,survive(s),&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1),&
                    shock_pi(sim_shocks(s,it),:,s,sim_types(s,it)),&
                    v(sim_types(s,it),:,:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),&
                    sim_labor_alloc(s,it-start_row+1),partprod,hbar)
            if (s==1) then
                sim_v_alloc(it-start_row+1)=dummy
            end if
            if (sim_labor_alloc(s,it-start_row+1)<hbar) then
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s,it-start_row+1),sim_capital_alloc(s+1,it-start_row+1),tr,rate_use,&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod)
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/pe
                sim_earnings_alloc(s,it-start_row+1)=w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod*sim_labor_alloc(s,it-start_row+1)
            else
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s,it-start_row+1),sim_capital_alloc(s+1,it-start_row+1),tr,rate_use,&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1))
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/pe
                sim_earnings_alloc(s,it-start_row+1)=w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*sim_labor_alloc(s,it-start_row+1)
            end if
                sim_labor_tax_alloc(s,it-start_row+1)=all_taxer(0.0D0,sim_earnings_alloc(s,it-start_row+1))
                sim_capital_tax_alloc(s,it-start_row+1)=lambdak*rate_use*(tr+sim_capital_alloc(s,it-start_row+1))
        end do


        do s=Jnr,J
            con=-1.0D0
            sim_labor_alloc(s,it-start_row+1)=0.0D0
            if(sim_shocks(Jnr,it) .le. shocksN/2) then
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1),&
                    kgrid(capital_positive:),&
                    tr,r,ss(sim_types(s,it),sim_shocks(s,it)),avg_earnings)
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))
                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it),sim_shocks(s,it)),tr,r,survive(s),&
                    v_retired(sim_types(s,it),sim_shocks(s,it),:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),avg_earnings,s)
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s+1,it-start_row+1),tr,r,ss(sim_types(s,it),sim_shocks(s,it)),avg_earnings)
            else
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1),&
                    kgrid(capital_positive:),&
                    tr,r,ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),avg_earnings)
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))
                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),tr,r,survive(s),&
                    v_retired(sim_types(s,it),sim_shocks(s,it),:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),avg_earnings,s)
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s+1,it-start_row+1),tr,r,ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),avg_earnings)
            end if
            dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
            sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/pe
            sim_capital_tax_alloc(s,it-start_row+1)=lambdak*r*(tr+sim_capital_alloc(s,it-start_row+1))

        end do
    end do

    call MPI_SEND( start_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( end_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_labor_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_earnings_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_capital_alloc,(J+1)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_consumption_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_energy_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_labor_tax_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_capital_tax_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_v_alloc,(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)


    deallocate(sim_labor_alloc)
    deallocate(sim_capital_alloc)
    deallocate(sim_energy_alloc)
    deallocate(sim_earnings_alloc)
    deallocate(sim_consumption_alloc)
    deallocate(sim_labor_tax_alloc)
    deallocate(sim_capital_tax_alloc)
    deallocate(sim_v_alloc)

    CALL MPI_Bcast( tax_diff, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    end do
call MPI_Bcast( diff_parameter,1,MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
end do




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Slaves that are not going to work for some
else if (my_id > (num_procs_run-1)) then
DO WHILE (diff_parameter>tol)
tax_diff=1.0D0
    DO WHILE (tax_diff >tol)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Solving Value Function (These are not working)

            CALL MPI_Bcast( prices, 9, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            CALL MPI_Bcast( ss, typesN*(shocksN/2), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            tr=prices(1)
            lambda0=prices(2)
            w=prices(3)
            r=prices(4)
            tauss=prices(5)
            ss_cap=prices(7)
            avg_earnings=prices(8)
            pe=prices(9)
            avg_earnings_scale=avg_earnings*(1.0D0-tauss/2.0D0)


        do s=Jnr-1,2,-1
            CALL MPI_Bcast( v_in(:,:,:), typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
       end do




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Simulating
    CALL MPI_Bcast( v(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    CALL MPI_Bcast( v_retired(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    call MPI_RECV(start_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
    call MPI_RECV(end_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

    allocate(sim_labor_alloc(J,end_row-start_row+1))
    allocate(sim_earnings_alloc(J,end_row-start_row+1))
    allocate(sim_capital_alloc(J+1,end_row-start_row+1))
    allocate(sim_consumption_alloc(J,end_row-start_row+1))
    allocate(sim_energy_alloc(J,end_row-start_row+1))
    allocate(sim_capital_tax_alloc(J,end_row-start_row+1))
    allocate(sim_labor_tax_alloc(J,end_row-start_row+1))
    allocate(sim_v_alloc(end_row-start_row+1))



    !!Initialize

    sim_capital_alloc(1,:)=0.0D0
    sim_labor_alloc=0.0D0
    sim_earnings_alloc=0.0D0
    sim_capital_tax_alloc=0.0D0
    sim_labor_tax_alloc=0.0D0





    do it=start_row,end_row
        !Working
        do s=1,Jnr-1
            if(sim_capital_alloc(s,it-start_row+1)<0.0D0) then
                if(s==1) then
                    rate_use=r/survive(1)
                else
                    rate_use=r/survive(s-1)
                end if
            else
                rate_use=r
            end if

            CALL working_valuer(s,sim_capital_alloc(s,it-start_row+1),kgridN,tr,rate_use,survive(s),&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1),&
                    shock_pi(sim_shocks(s,it),:,s,sim_types(s,it)),&
                    v(sim_types(s,it),:,:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),&
                    sim_labor_alloc(s,it-start_row+1),partprod,hbar)
            if (s==1) then
                sim_v_alloc(it-start_row+1)=dummy
            end if
            if (sim_labor_alloc(s,it-start_row+1)<hbar) then
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s,it-start_row+1),sim_capital_alloc(s+1,it-start_row+1),tr,rate_use,&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod)
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/pe
                sim_earnings_alloc(s,it-start_row+1)=w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod*sim_labor_alloc(s,it-start_row+1)
            else
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s,it-start_row+1),sim_capital_alloc(s+1,it-start_row+1),tr,rate_use,&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1))
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/pe
                sim_earnings_alloc(s,it-start_row+1)=w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*sim_labor_alloc(s,it-start_row+1)
            end if
                sim_labor_tax_alloc(s,it-start_row+1)=all_taxer(0.0D0,sim_earnings_alloc(s,it-start_row+1))
                sim_capital_tax_alloc(s,it-start_row+1)=lambdak*rate_use*(tr+sim_capital_alloc(s,it-start_row+1))
        end do


        do s=Jnr,J
            sim_labor_alloc(s,it-start_row+1)=0.0D0
            con=-1.0D0
            if(sim_shocks(Jnr,it) .le. shocksN/2) then
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1),&
                    kgrid(capital_positive:),&
                    tr,r,ss(sim_types(s,it),sim_shocks(s,it)),avg_earnings)
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))

                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it),sim_shocks(s,it)),tr,r,survive(s),&
                    v_retired(sim_types(s,it),sim_shocks(s,it),:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),avg_earnings,s)
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s+1,it-start_row+1),tr,r,ss(sim_types(s,it),sim_shocks(s,it)),avg_earnings)
            else
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1),&
                    kgrid(capital_positive:),&
                    tr,r,ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),avg_earnings)
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))

                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),tr,r,survive(s),&
                    v_retired(sim_types(s,it),sim_shocks(s,it),:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),avg_earnings,s)
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s+1,it-start_row+1),tr,r,ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),avg_earnings)
            end if
            dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
            sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/pe
            sim_capital_tax_alloc(s,it-start_row+1)=lambdak*r*(tr+sim_capital_alloc(s,it-start_row+1))
!            sim_capital_tax_alloc(s,it-start_row+1)=ret_taxer(sim_capital_alloc(s,it-start_row+1),&
!                sim_capital_alloc(s+1,it-start_row+1),tr,r,ss,avg_earnings)
        end do
    end do

    call MPI_SEND( start_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( end_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_labor_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_earnings_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_capital_alloc,(J+1)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_consumption_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_energy_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_labor_tax_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_capital_tax_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_v_alloc,(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)


    deallocate(sim_labor_alloc)
    deallocate(sim_capital_alloc)
    deallocate(sim_energy_alloc)
    deallocate(sim_earnings_alloc)
    deallocate(sim_consumption_alloc)
    deallocate(sim_labor_tax_alloc)
    deallocate(sim_capital_tax_alloc)
    deallocate(sim_v_alloc)

    CALL MPI_Bcast( tax_diff, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    end do
call MPI_Bcast( diff_parameter,1,MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
end do





end if
deallocate(population)








call MPI_FINALIZE(ierr)

end PROGRAM hetero_ss
